#ifdef _DEBUG

#include <cassert>

using namespace matrix;

void testVect_prod() {
    // prepare
    TMatrix<double, 3, 1> a;
    TMatrix<double, 1, 3> b;
    a(0) = 1;    b(0, 0) = 4;   b(0, 1) = 5;   b(0, 2) = 6;
    a(1) = 2;
    a(2) = 3;
    // act
    TMatrix<double, 3, 3> c = a * b;
    // check
    assert(c(0, 0) = 4);    assert(c(1, 0) = 8);    assert(c(2, 0) = 12);
    assert(c(0, 1) = 5);    assert(c(1, 1) = 10);   assert(c(2, 1) = 15);
    assert(c(0, 2) = 6);    assert(c(1, 2) = 12);   assert(c(2, 2) = 18);
    // prepare
    TMatrix<double, 3, 1> x;
    TMatrix<double, 3, 1> y;
    x(0) = 1;   y(0) = 4;
    x(1) = 2;   y(1) = 5;
    x(2) = 3;   y(2) = 6;
    // act
    double z = prod(x, y);
    // check
    assert(z == 4 + 10 + 18);
}

void testTrans() {
    // prepare
    TMatrix<double, 3, 5> a;
    size_t m, n;
    for(m = 0; m < a.size1; m++) {
        for(n = 0; n < a.size2; n++) {
            a(m, n) = m * 10 + n;
        }
    }
    // act
    const TMatrix<double, 5, 3> b = trans(a);
    // check
    for(m = 0; m < b.size1; m++) {
        for(n = 0; n < b.size2; n++) {
            assert(b(m, n) == n * 10 + m);
        }
    }
}

void testAdd() {
    // prepare
    TMatrix<double, 3, 5> a;
    TMatrix<double, 3, 5> b;
    size_t m, n;
    for(m = 0; m < a.size1; m++) {
        for(n = 0; n < a.size2; n++) {
            a(m, n) = m * 10. + 1.*n;
            b(m, n) = m * 1000. + n * 100.;
        }
    }
    // act
    const TMatrix<double, 3, 5> c = b + a;
    const TMatrix<double, 3, 5> d = b - a;
    // check
    for(m = 0; m < c.size1; m++) {
        for(n = 0; n < c.size2; n++) {
            assert(c(m, n) == (m * 1000 + n * 100) + (m * 10 + n));
            assert(d(m, n) == (m * 1000 + n * 100) - (m * 10 + n));
        }
    }
}

void testLong_expressions() {
    // prepare
    TMatrix<double, 3, 1> a;
    TMatrix<double, 1, 3> b;
    TMatrix<double, 3, 3> c;
    TMatrix<double, 3, 3> d;
    a(0) = 1;    b(0, 0) = 4;
    a(1) = 2;    b(0, 1) = 5;
    a(2) = 3;    b(0, 2) = 6;
    for(size_t m = 0; m < c.size1; m++)
        for(size_t n = 0; n < c.size2; n++)
            c(m, n) = m * n + n;
    // act
    d = c + a * b;
    // check
    assert(d(0, 0) == 4);    assert(d(1, 0) == 8);    assert(d(2, 0) == 12);
    assert(d(0, 1) == 6);    assert(d(1, 1) == 12);   assert(d(2, 1) == 18);
    assert(d(0, 2) == 8);    assert(d(1, 2) == 16);   assert(d(2, 2) == 24);
}

void testTriangular() {
    // prepare
    TMatrix<double, 3, 3> d;
    d(0, 0) = 4;    d(1, 0) = 8;    d(2, 0) = 12;
    d(0, 1) = 6;    d(1, 1) = 12;   d(2, 1) = 18;
    d(0, 2) = 8;    d(1, 2) = 16;   d(2, 2) = 24;
    // act
    TTriangular<double, 3, 3> tri(d);
    // check
    assert(tri(0, 0) == 4);    assert(tri(1, 0) ==  8);   assert(tri(2, 0) == 12);
    assert(tri(0, 1) == 0);    assert(tri(1, 1) == 12);   assert(tri(2, 1) == 18);
    assert(tri(0, 2) == 0);    assert(tri(1, 2) ==  0);   assert(tri(2, 2) == 24);
}

void testCholesky_invert() {
    // prepare
    TMatrix<double, 3, 3> d;
    d(0, 0) = 1;    d(0, 1) = 1;   d(0, 2) = 1;
    d(1, 0) = 1;    d(1, 1) = 1;   d(1, 2) = 1;
    d(2, 0) = 1;    d(2, 1) = 1;   d(2, 2) = 1;
    // act
    TMatrix<double, 3, 3> i = cholesky_invert(d);
    // check
    assert(i(0, 0) == 2);    assert(i(1, 0) ==-1);   assert(i(2, 0) == 0);
    assert(i(0, 1) ==-1);    assert(i(1, 1) == 2);   assert(i(2, 1) ==-1);
    assert(i(0, 2) == 0);    assert(i(1, 2) ==-1);   assert(i(2, 2) == 1);
}

void testColumn_and_row() {
    // prepare
    TMatrix<double, 3, 3> a;
    a(0, 0) = 1;    a(0, 1) = 2;   a(0, 2) = 3;
    a(1, 0) = 1;    a(1, 1) = 2;   a(1, 2) = 3;
    a(2, 0) = 1;    a(2, 1) = 2;   a(2, 2) = 3;
    // act
    TMatrix<double, 3, 1> c1 = column(a, 0);
    TMatrix<double, 3, 1> c2 = column(a, 2);
    TMatrix<double, 3, 1> d1 = c1 + c2;
    // check
    assert(c1(0) == 1);    assert(c1(1) == 1);   assert(c1(2) == 1);
    assert(c2(0) == 3);    assert(c2(1) == 3);   assert(c2(2) == 3);
    assert(d1(0) == 4);    assert(d1(1) == 4);   assert(d1(2) == 4);
    // act
    TMatrix<double, 3, 1> r1 = row(a, 0);
    TMatrix<double, 3, 1> r2 = row(a, 2);
    TMatrix<double, 3, 1> d2 = r1 + r2;
    // check
    assert(r1(0) == 1);    assert(r1(1) == 2);   assert(r1(2) == 3);
    assert(r2(0) == 1);    assert(r2(1) == 2);   assert(r2(2) == 3);
    assert(d2(0) == 2);    assert(d2(1) == 4);   assert(d2(2) == 6);
}

void testCholesky_decompose() {
    // prepare
    TMatrix<double, 3, 3> a;
    a(0, 0) = 1;    a(0, 1) = 2;   a(0, 2) = 3;
    a(1, 0) = 2;    a(1, 1) = 8;   a(1, 2) = 2;
    a(2, 0) = 3;    a(2, 1) = 2;   a(2, 2) = 30;
    // act
    TMatrix<double, 3, 3>l = cholesky_decompose(a);
    // check
    TMatrix<double, 3, 3>b = prod(l, trans(l));
    for(size_t i = 0; i < a.size1; i++) {
        for(size_t j = 0; j < a.size2; j++) {
            assert(a(i, j) == b(i, j));
        }
    }
}

void testInvert() {
    // Testing fact:
    //    A * (A**-1) == I
    // prepare
    TMatrix<double, 3, 3> a;
    a(0, 0) = 1;    a(0, 1) = 2;   a(0, 2) = 3;
    a(1, 0) = 2;    a(1, 1) = 8;   a(1, 2) = 2;
    a(2, 0) = 3;    a(2, 1) = 2;   a(2, 2) = 30;
    // act
    TMatrix<double, 3, 3>inv = invert(a);
    // check
    TMatrix<double, 3, 3>I = prod(a, inv);
    TMatrix<double, 3, 3>ident;
    ident.identity();
    double epsilon = 1.e-15;
    for(size_t i = 0; i < a.size1; i++)
        for(size_t j = 0; j < a.size2; j++)
            assert( abs(I(i, j) - ident(i, j)) < epsilon );
}

namespace matrix {
    void testAll() {
        testTriangular();
        testColumn_and_row();
        testTrans();
        testAdd();
        testVect_prod();
        testLong_expressions();
        testCholesky_decompose();
        testCholesky_invert();
        testInvert();
    }
}

#endif